Doublon relaxation in the Bose-Hubbard model 
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Decay of a high-energy double occupancy state, doublon, in a narrow-band lattice requires creation 
of a coherent many-particle excitation. This leads to an exponentially long relaxation time of such 
a state. We show that, if the average occupation number is sufficiently small, the corresponding 
exponent may be evaluated exactly. To this end we develop the quasiclassical approach to calculation 
of the high-order tree-level decay amplitudes. 



Recent experiments with cold atomic gases in opti- 
cal lattices [1] allowed to create realizations of Hubbard 
model with high degree of control over the main parame- 
ters, the bandwidth W and the on-site interaction U. In 
addition to exploring rich equilibrium phase diagram [5] , 
a rapid time variation of the parameters allows to prepare 
atomic systems in highly excited states suitable for stud- 
ies of the far from equilibrium dynamics of this strongly 
correlated system. The most prominent high-energy exci- 
tations are repulsively bound doubly occupied sites called 
doublons [5] . They were recently observed in experiments 
with both bosonic [5] and fcrmionic [I] atoms. The dy- 
namical properties of doublons were considered in recent 
publications [5H2]- 

In experiments doublon energy Ud = U + W stays be- 
tween the first and second Bloch bands, which allows to 
consider the single band Hubbard model. In the absence 
of other particles doublon is the stable excitation with 
a very heavy mass pJj- If other particles with the aver- 
age occupation po are present, they offer a possibility for 
doublon to decay by transferring its interaction energy 
to the kinetic energy of single-particle excitations. In the 
interesting limit v = Ud/W 1, the number of such final 
state excited particles n > v is large. Calculation of the 
relaxation rate requires therefore analyzing very high or- 
ders of the perturbation theory in the inter-particle inter- 
actions. Such an analysis for nearly half-filled fcrmionic 
model was recently performed in Ref. [5], which found 
that the decay rate scales as r _1 oc exp{— av In v}, where 
a was found to be approximately 0.8. 

In this letter we show that if the problem admits a 
small parameter - the average filling factor po <C 1, 
the leading exponential scaling of the decay rate may be 
found exactly. We use the Bose-Hubbard version of the 
model for the illustration. Our approach allows one to 
calculate generating function of all n-particle tree-level 
threshold (i.e. such that v = n) amplitudes. Such a gen- 
crating function is shown to obey classical equation of 
motion for the coupled doublon and particle fields. Luck- 
ily this equation admits an analytic solution, which yields 
exact expression for the tree-level amplitudes. The use of 
the tree approximation is justified as long as p§v < 7^, 



where 7d = e~ 3 (Tr/2e) D / 2 and D is spatial dimensional- 
ity of the lattice. The similar approach was developed in 
high-energy physics for calculation of the threshold am- 
plitudes for 1 — > n and 2 — > n deep inelastic scattering 
processes [51 H]. 

If po <C 7_d, the window of parameters 1 <C v < Jd/Po 
exists where the loop diagrams may be neglected as 
carrying extra small factors po and only tree diagrams 
contribute to the decay rate. They may be effectively 
summed up with the quasiclassical procedure, described 
below, leading to the doublon decay rate 
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The pre-exponential factor is affected by the loop dia- 
grams and is beyond the accuracy of the quasiclassical 
calculation. The factor [(2/D)ln(^/r)/p v)] D / 2 originates 
from spreading wavepacket of excited particles in D spa- 
tial dimensions [10] . Apart from this factor, the result in 
Eq. ([I]) reminds Landau expression for the multi-quanta 
excitation of a nonlinear zero-dimensional oscillator |11) . 

We consider Bose-Hubbard model on a D-dimensional 
cubic lattice 
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where (i, j) denote nearest neighbor sites. The first term 
in the Hamiltonian, Hq describes the hopping between 
adjacent sites and creates a band of the delocalized sin- 
gle particle states of the width W — AJD proportional to 
the tunneling amplitude J pj [2] . Its dispersion relation 
is symmetric around zero energy ep, = — 2J^u=i cosfc M , 
where — n < < n is the quasimomentum within the 
first Brillouin zone, measured in units of the inverse lat- 
tice spacing a. In the optical lattice, the lattice spacing 
is given by the half of the laser wave length, while the 
on-site repulsion U is proportional to the bulk s-wave 
scattering length a s [HE], which in turn can be tuned 
experimentally in the wide range by taking advantage of 
Feshbach resonance phenomenon. 
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The zero-temperature phase diagram of this model 
consists of the superfiuid and the Mott-insulator phases 
[IJ [12] . The Mott insulator is characterized by an integer 
occupation number p at each site and a finite energy 
gap for the quasiparticle excitations. On the other hand, 
the compressible superfiuid phase is characterized by the 
presence of the Bose condensate with a noninteger aver- 
age occupation p$. We restrict ourselves with the latter 
scenario at po <C 1. 

Among the momentum content Ck of the bosonic fields 
a special role is played by the field co(t) = y/p~oe~' lflt , 
representing the condensate with the chemical potential 
p = —W/2. We shall also separate the field c v (t), de- 
scribing quantized particles with the energy e n — p = W 
at the very top of the Bloch band. To proceed we pass 
to a coherent state functional representation for the evo- 
lution operator e lHt and decouple the interation term in 
Eq. ([2| with the help of the Hubbard-Stratonovich trans- 
formation 

2 llll ~2XJ 2 [l,l \d l Oj\ci)' 1 ' 

introducing the auxiliary field di(t) which describes the 
doublon. We now perform Gaussian integration over 
bosonic fields c k with k ^ 0, 7r, using the free propagator 
Gfc(e) = (e — and expand the resulting action to 

the second order in the d- fields [13]. This leads to the 
doublon Lagrangian 

^ = \ E [ U ~ l - C ^ <*)K( e ) , (4) 
q 

where C(e,q) is the Cooper polarization given by 

C(e, q) = i£ G k+q {e - e')G. k {e') = £ ^ • 

(5) 

The poles of the doublon propagator, i.e. C(e,q) = U 1 , 
determine its dispersion relation e = e^{q). In the limit 
U J one may expand the denominator on the right 
hand side of Eq. ^ to the second order to find [3] 

oj2 D „ 

Q(g)^ + ^£cos 2 | + o(j 4 /c/ 3 ). (6) 

Therefore the doublon band is narrow W 2 /2DU <C W. 
One can thus disregard the q-dependence of the doublon 
dispersion and think of it as of infinitely heavy localized 
particle with the energy U. This is equivalent to the 
approximation where one takes C(e,q) ps C(e, n) = er 1 . 
Hereafter we suppress site or momentum index of the 
doublon, assuming it to be localized and infinitely heavy. 

In the absence of the condensate, po = 0, the doublon 
is absolutely stable, which is reflected in the presence of 
the true pole in its Green function, cf. Eq EL at e ss U. 
Interaction with the condensate results in the doublon 




FIG. 1. An example of a tree- level diagram for the ampli- 
tude of the doublon decay in n = 5 particles, a) Doublon 
propagator (dd); b) particle propagator (even-); c) condensate 
amplitude Co oc ^/po\ d) doublon decay vertex c^-c^-d; e) par- 
ticle interaction with the condensate vertex dcoc n . Crossed 
lines denote the mass-shell particles. 



self-energy E(e). Its imaginary part taken at the mass- 
shell is the doublon decay rate l/2r = ImS(J7), which is 
the focus of this letter. With the help of the vertex c^cvd, 
cf. Eq. pi), the doublon initially decays onto two virtual 
particles with momenta close to the boundaries of the 
Brillouin zone. Each of the created particles may collide 
with the condensate and with the help of the vertex dc v Co 
create more virtual doublons, the latter again decay, etc. 
This process, depicted in Fig.JT] continues until the initial 
doublon energy (counted from the chemical potential) 
Ud = U — 2p = U + W is transferred into the kinetic 
energy of n > Ud/W real particles. The corresponding 
decay rate is given by the Golden Rule 

\ = E Y E \A n \ 2 s( Ud -±(e pi -p)\, (7) 

n>U d /W Pl,-,Pn \ 1=1 ) 



An(U d ;pi,.-.Pn) = (pi,...,Pn|44d|d), 

where (p% ,.. .p n \ is the final n-particle excited state of 
the interacting gas and \d) is the initial state with the 
single doublon and rest of the gas in the condensate. 

Consider now a special situation where the doublon 
energy happens to be slightly below the n-particle cre- 
ation threshold Ud < nW . In this case all n final 
particles must be very close to the top of the band. 
One may thus approximate their dispersion relation as 
e Pl — p = W — (pi — it) 2 /2m*, where m* = (2J)^ 1 is the 
effective mass. The corresponding phase volume for the 
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decay process 



1=1 



[(nW-U d )/4TrJ}^ 
T(^)(nW-U d ) 



(8) 

depends very strongly on the deviation nW — Ud > 
from the threshold. On the other hand, the correspond- 
ing matrix element A n {Ud',Pi, ■ ■ -Pn) appears to be a 
smooth function of Ud- This implies that in the vicin- 
ity of the threshold it may be substituted by a constant 
A n ~ A n (nW; ir, . . .71"). Indeed, once the doublon en- 
ergy is fixed at the n-particle threshold Ud = nW, all 
final state particles must be at the top of the band, i.e. 
in the same quantum state with p — n. Therefore only 
p = 7r state (along with the condensate p = 0) may ever 
appear as a real particle, while the rest of (already in- 
tegrated out) states are purely virtual. As a result the 
decay rate of the doublon with the energy Ud is given by 

! ^ E K| 2 V„(£/ d ). (9) 

n>U d /W 
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The crack of the matter thus is to evaluate the threshold 
amplitudes A n . 

If the doublon energy Ud is not too large the leading 
contribution to the amplitude comes from the tree-level 
diagrams, i.e. those which do not include loops. Indeed 
for the same number of the final state mass-shell particles 
n, a loop diagram involves a higher order of the perturba- 
tion theory than a corresponding tree one. Such tree level 
diagrams are generated by repeated Wick contractions of 
the two types of vertices c^c^d and dc^c^. The Hermitian 
conjugated vertices participate only in the loop diagrams. 
As a result, the tree level amplitudes are generated by the 
following Lagrangian 
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(10) 

where the condensate field possesses a classical expecta- 



tion value co(t) — J~p§e ' . Since the tree-generating 
Lagrangian is linear in annihilation fields c n (t) and d(t), 
the functional integration over these fields enforces equa- 
tions of motion for the corresponding creation compo- 
nents c n (t) and d(t). To avoid the time non-local op- 
erator 1/e let us define an auxiliary dimensionless field 
a = d/e, i.e. d(t) = —idt<j{t). It is also convenient to 
shift the energies for the condensate to be at zero by the 
gauge transformation c n — > c^e~ lWt ^ 2 and a — > ae~ 



-iWt 



Taking then variation of the Lagrangian ( 10 1 with re- 



spect to Ctt and d, one obtains the following equations of 
motion 

(—id t — W) = y/p~o {—idt — W) a ; (11) 
{-id t -U d )a = Uc n c w . (12) 

The solution we seek may be written as a series of positive 
powers of z{t) = e lWt . For example c % {t) = Y^=i a nZ n , 



where a n represent the tree diagrams starting with one 
virtual particle and ending up with n mass-shell particles, 
each having the energy W (thus the factor e mWt ). The 
proper normalization is thus cx\ = 1. Similarly a(t) = 
Y^n°=i PnZ n , generates tree diagrams which start from the 
doublon and end up with n mass-shell particles, Fig. [T] 
Since the doublon must decay at least on two particles, 
the normalization is (3\ = 0. 



The properly normalized solution of Eq. ( 11 ) is thus 
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Expressing a and employing Eq. (12), one finds 
(id t + U d )c w + U^cl = Ue mt . 



(13) 



(14) 



This equation describes a weakly non-linear high fre- 
quency Ud oscillator, which is forced with the low fre- 
quency W U d driving force. The non-linearity gener- 
ates higher harmonics of the applied force, bringing the 
oscillator into the exact resonance if U d /W = n is an 
integer. One expects thus the solution to have a pole 
(nW — Ud)" 1 , reflecting the resonance condition. This 
pole represents the Green function of the incoming dou- 
blon at the threshold energy e = nW . Since we are in- 
terested in the self-energy Y,(nW) rather than the Green 
function itself, we need to focus only on the residue of 
the corresponding pole. The latter is a function of time 
given by a series in powers of z. The threshold am- 
plitude A n is proportional to the coefficient in front of 
term of this series. Indeed, it is exactly the 
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e inW{t-t) in (d(t)d(t')) , which generates 
the proper delta-function in Eq. ([T]) upon the Fourier 
transform. 



Equation (14) is of Riccati type, which may be trans- 



formed into the linear second order differential equation 
by the substitution = idtv / (U ^/po v) . Changing also 
the variable t — > z and using v = Ud/W, one finds 
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0. 



(15) 



The exact solution of Eq. (15) is given in terms of the 
Bessel functions 

v(z) = (bz) v ' 2 [CxT{l - v)I- v {2y/bz)+ 

(-1)^2^1 + ^(2^)}, (16) 

where b — ^fp${y — l) 2 , and C\ and C2 are free inte- 
gration constants. As explained above, we look for the 
resonant poles of the form z n j(n — v) at integer values 
of v. It is easy to see that they may come only from the 
first term on the r.h.s. of Eq. (Tl6|), which reads as 



v{z) = C 1 T{1-v)Y j 



(bz) h 
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The pole structure can be made explicit by rewriting the 
coefficient of z k in the form 



r(i-i,) 



t 



r(jfe + 1-1/) (k - v)(k - 1 - v){k - 2 - v) ... (1 - v) 

(18) 

At integer v = n it contains poles for k = n, n + 1, . . .. 
The resonant time dependence z n is provided only by 
the pole with the lowest value k = n. Other pole terms 
contribute to the decay processes with more than n out- 
going particles and are neglected. Notice also that at 
z — >• one has = —(W/Uy/po)zd z bxv = z + 0(z 2 ), 
which provides the proper normalization of the generat- 
ing function for tree diagrams. The appropriate term in 
the expansion of c^(z) takes the form [c K ] n z n / (n — u), 
where the coefficient 
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is ev alua ted at v = n, Recalling that &(z) — (cv— z)/ yfpo, 
Eq. (131, while the doublon field is d(z) = Wzd z a(z), 
one finds for the n-particle threshold amplitude A n = 



An = WVn\ (-1) 



„ ( n - 1) 2 >| 
[(n-1)!] 2 



(20) 



responding to f = 1. This amounts to i w 2.25. Em- 
ploying the same functional form of the prefactor as in 
Eq. Q we estimate fa. (2.25) « 4.5. 

In conclusion, the decay of doublon is accompanied 
by creation of a many-particle excitation. In a certain 
range of parameters, dictated by a small filling factor, 
such process is described by tree-level diagrams. If this 
is the case, it is amenable to the semiclassical evaluation 
close in spirit to Ref. [11] . We performed this calcula- 
tion for the Bose-Hubbard model in the superfluid phase 
and determined exactly the leading terms in the decay 
rate exponent. It remains to be seen if the method and 
the results can be adopted to the Fermi-Hubbard model, 
treated in Ref. [5] by other means. 
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where the proportionality coefficient yn\ originates from 
the n! ways of pairing n final on mass-shell cv-particlcs 
in the doublon self-energy. 

The total decay rate of the doublon is calculated by 
summing up over all open decay channels according to 
Eq. ([9]). Substituting explicit expressions for the decay 
amplitude Eq. (20) and the phase volume Eq. pi), we 



obtain the decay rate in form of the asymptotic series 
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To sum up the series one employs Stirling formula for 
gamma functions, substitutes summation by the integra- 
tion and performs the latter in the stationary point ap- 
proximation. As a result, one obtains the expression ([!]) 
for the doublon decay rate. One may argue that the 
most general scaling form of the decay time is given by 
lnr = pQ 1 hu(x), where x = p^v. The dimensionless 
function ho{x) is known in high energy physics as a 
"Holy Grail" function [5] . The semiclassical methods 
allow to evaluate it in the limit x <C 1, cf. Eq. ([lj. There 
are arguments (see Ref. [8]) that it either saturates, or 
grows extremely slowly at x > 1. Its behavior may be 
estimated from comparison with experimental results re- 
ported in Ref. [3]. The reported value of relaxation time 
t s» 700ms was measured for v w 7.5, po ~ 0.3, and 
W/h rj 1583Hz in highly anisotropic optical lattice cor- 
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